
R version 4.5.2 (2025-10-31) -- "[Not] Part in a Rumble"
Copyright (C) 2025 The R Foundation for Statistical Computing
Platform: aarch64-apple-darwin20

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> # Marijuana experiment
> 
> # Load relevant libraries
> library(data.table)
> library(tidyverse)
-- Attaching core tidyverse packages ------------------------ tidyverse 2.0.0 --
v dplyr     1.1.4     v readr     2.1.6
v forcats   1.0.1     v stringr   1.6.0
v ggplot2   4.0.1     v tibble    3.3.0
v lubridate 1.9.4     v tidyr     1.3.1
v purrr     1.2.0     
-- Conflicts ------------------------------------------ tidyverse_conflicts() --
x dplyr::between()     masks data.table::between()
x dplyr::filter()      masks stats::filter()
x dplyr::first()       masks data.table::first()
x lubridate::hour()    masks data.table::hour()
x lubridate::isoweek() masks data.table::isoweek()
x dplyr::lag()         masks stats::lag()
x dplyr::last()        masks data.table::last()
x lubridate::mday()    masks data.table::mday()
x lubridate::minute()  masks data.table::minute()
x lubridate::month()   masks data.table::month()
x lubridate::quarter() masks data.table::quarter()
x lubridate::second()  masks data.table::second()
x purrr::transpose()   masks data.table::transpose()
x lubridate::wday()    masks data.table::wday()
x lubridate::week()    masks data.table::week()
x lubridate::yday()    masks data.table::yday()
x lubridate::year()    masks data.table::year()
i Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
> library(broom)
> library(Cairo)
> library(ggthemes)
> library(knitr)
> library(kableExtra)

Attaching package: 'kableExtra'

The following object is masked from 'package:dplyr':

    group_rows

> library(psych)

Attaching package: 'psych'

The following objects are masked from 'package:ggplot2':

    %+%, alpha

> library(xtable)
> 
> red_mit <- '#A31F34'
> red_light <- '#A9606C'
> blue_mit <- '#315485'
> 
> ############################
> # Data Processing
> ############################
> 
> w1 <- read_csv("Data/MediaSSI_Dec2017_w1_recoded.csv") %>%
+   filter(consent == 1)
Rows: 7394 Columns: 35
-- Column specification --------------------------------------------------------
Delimiter: ","
chr  (7): StartDate, EndDate, ResponseId, med_pref, med_choice, PID, article...
dbl (27): X, Progress, consent, gender, educ, income, party1, party4, forced...
num  (1): race

i Use `spec()` to retrieve the full column specification for this data.
i Specify the column types or set `show_col_types = FALSE` to quiet this message.
> w2 <- read_csv("Data/MediaSSI_Dec2017_w2_recoded.csv")
Rows: 4927 Columns: 24
-- Column specification --------------------------------------------------------
Delimiter: ","
chr  (7): StartDate, EndDate, ResponseId, med_pref, med_choice, PID, article...
dbl (17): Progress, msnbc, fox, entertainment, forcedchoice, pid, mar_costmo...

i Use `spec()` to retrieve the full column specification for this data.
i Specify the column types or set `show_col_types = FALSE` to quiet this message.
> w3 <- read_csv("Data/MediaSSI_Dec2017_w3_recoded.csv")
Rows: 4527 Columns: 17
-- Column specification --------------------------------------------------------
Delimiter: ","
chr  (5): StartDate, EndDate, ResponseId, med_pref, PID
dbl (12): Progress, mar_costmore, mar_fewserious, mar_wrong, mar_violence, m...

i Use `spec()` to retrieve the full column specification for this data.
i Specify the column types or set `show_col_types = FALSE` to quiet this message.
> 
> keep_best <- \(df, prog) df %>%
+   filter(!is.na(.data[[prog]])) %>%
+   arrange(PID, desc(.data[[prog]])) %>%
+   group_by(PID) %>% slice_head(n = 1) %>% ungroup()
> 
> w1u <- keep_best(w1, "Progress")
> w2u <- keep_best(w2, "Progress")
> w3u <- keep_best(w3, "Progress")
> 
> w123 <- w1u %>%
+   left_join(w2u, by = "PID", suffix = c("", "_w2")) %>%
+   left_join(w3u, by = "PID", suffix = c("", "_w3"))
> 
> 
> 
> # Define the outcomes
> outcomes <- c("mar_tradeoff", "mar_econ", "mar_costmore", "mar_fewserious",
+               "mar_wrong", "mar_violence", "mar_legmed", "mar_serious",
+               "mar_legrec", "danger_mar")
> 
> # Calculate Cronbach's alpha
> index_fa <- psych::alpha(w123[, outcomes], check.keys = TRUE)
> alpha <- index_fa$total["raw_alpha"]
> print(paste("Cronbach's alpha:", round(alpha, 2)))
[1] "Cronbach's alpha: 0.89"
> 
> w123 <- w123 %>%
+   mutate(
+     mar_index = rowMeans(.[outcomes], na.rm = TRUE),
+     mar_index_w2 = rowMeans(.[paste0(outcomes, "_w2")], na.rm = TRUE),
+     mar_index_w3 = rowMeans(.[paste0(outcomes, "_w3")], na.rm = TRUE)
+   )
> 
> calculate_pca <- function(data, vars) {
+   result <- rep(NA, nrow(data))
+   complete_cases <- complete.cases(data[vars])
+   if(sum(complete_cases) > 0) {
+     complete_data <- data[complete_cases, vars]
+     pca <- prcomp(complete_data, center = TRUE, scale. = TRUE)
+     pca_scores <- pca$x[,1]
+     
+     # Check correlation with the additive index
+     additive_index <- rowMeans(complete_data, na.rm = TRUE)
+     if(cor(pca_scores, additive_index) < 0) {
+       pca_scores <- -pca_scores  # Flip the sign if negatively correlated
+     }
+     
+     pca_scores <- scale(pca_scores)  # Standardize
+     result[complete_cases] <- pca_scores
+   }
+   return(result)
+ }
> 
> w123$pcw1 <- calculate_pca(w123, outcomes)
> w123$pcw2 <- calculate_pca(w123, paste0(outcomes, "_w2"))
> w123$pcw3 <- calculate_pca(w123, paste0(outcomes, "_w3"))
> 
> ###############################################
> # Figure 3 - Persuasive effects in full sample
> ###############################################
> 
> extract_coefs <- function(model) {
+   coefs <- coef(model)
+   se <- sqrt(diag(vcov(model)))
+   data.frame(
+     estimate = coefs[-1],
+     se = se[-1],
+     treatment = c("Fox", "MSNBC"),
+     conf.low = coefs[-1] - qnorm(0.975) * se[-1],
+     conf.high = coefs[-1] + qnorm(0.975) * se[-1],
+     conf.low.90 = coefs[-1] - qnorm(0.95) * se[-1],
+     conf.high.90 = coefs[-1] + qnorm(0.95) * se[-1]
+   )
+ }
> 
> # First Principal Component Plot
> pc_w1 <- lm(pcw1 ~ fox + msnbc, data = filter(w123, forcedchoice == 1 & forcedchoice_w2 == 0))
> pc_w2 <- lm(pcw2 ~ fox + msnbc, data = filter(w123, forcedchoice == 1 & forcedchoice_w2 == 0))
> pc_w3 <- lm(pcw3 ~ fox + msnbc, data = filter(w123, forcedchoice == 1 & forcedchoice_w2 == 0))
> 
> fig1_data_pc <- bind_rows(
+   extract_coefs(pc_w1) %>% mutate(wave = 1),
+   extract_coefs(pc_w2) %>% mutate(wave = 2),
+   extract_coefs(pc_w3) %>% mutate(wave = 3)
+ ) %>% mutate(outcome = "First Principal Component")
> 
> labs_pc <- data.frame(
+   wave=c(1.5, 1.5), 
+   estimate=c(0.25, -0.2), 
+   label=c("Fox rather\nthan Entertainment", "MSNBC rather\nthan Entertainment"),
+   treatment=c("Fox", "MSNBC")
+ )
> 
> pdf('Output/fig3_pc.pdf', width=6, height=4)
> print(ggplot(fig1_data_pc, aes(y=estimate, x=wave, color=treatment)) +
+   geom_hline(yintercept = 0, lty=2) + 
+   geom_point(position=position_dodge(width=0.25))+ 
+   geom_line(aes(group=treatment), position=position_dodge(width=0.25)) + 
+   geom_errorbar(aes(ymin=conf.low, ymax=conf.high, width=0), size=0.5, position=position_dodge(width=0.25)) + 
+   geom_errorbar(aes(ymin=conf.low.90, ymax=conf.high.90, width=0), size=1, position=position_dodge(width=0.25)) + 
+   geom_text(data=labs_pc, aes(x=wave, y=estimate, label=label, color=treatment)) +
+   theme_bw() + 
+   scale_x_continuous('Survey Wave', breaks=c(1,2,3), limits=c(0.75,3.25)) + 
+   scale_y_continuous('Effect of Treatment on\nFirst Principal Component',
+                      breaks=seq(-0.5,0.5,0.25),
+                      labels=c("\n\n-0.5\n(More\nliberal)","-0.25",  "0.0","0.25", "(More\nconservative)\n0.5\n\n"),
+                      limits=c(-0.55,0.55)) + 
+   scale_color_manual("Treatment", breaks = c("Fox","MSNBC"), values=c(red_mit,blue_mit)) +
+   theme(legend.position="none", strip.background = element_rect(fill="lightgrey"), axis.title.y = element_text(margin = margin(r = -20))))
Warning message:
Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
i Please use `linewidth` instead. 
> dev.off()
null device 
          1 
> 
> 
> ## Figure 3 notes
> 
> # Analysis sample
> ana <- w123 %>% dplyr::filter(forcedchoice == 1 & forcedchoice_w2 == 0)
> 
> ctrl_mean_sd_N <- function(y) {
+   z <- ana[[y]][ ana$entertainment == 1 ]
+   c(mean = mean(z, na.rm = TRUE), sd = sd(z, na.rm = TRUE), N = sum(!is.na(z)))
+ }
> 
> 
> # --- First PC: control-group means/SD (and Ns) by wave
> pc_w1_ctrl <- ctrl_mean_sd_N("pcw1")
> pc_w2_ctrl <- ctrl_mean_sd_N("pcw2")
> pc_w3_ctrl <- ctrl_mean_sd_N("pcw3")
> 
> # Formatters
> fmt  <- function(x) sprintf("%.2f", x)
> fmt0 <- function(x) format(round(x), big.mark = ",")
> 
> # Build caption strings
> 
> caption_pc <- paste0(
+   "Control (Entertainment) means [SD], N — Wave 1: ", fmt(pc_w1_ctrl["mean"]), " [", fmt(pc_w1_ctrl["sd"]), "], N=", fmt0(pc_w1_ctrl["N"]), "; ",
+   "Wave 2: ", fmt(pc_w2_ctrl["mean"]), " [", fmt(pc_w2_ctrl["sd"]), "], N=", fmt0(pc_w2_ctrl["N"]), "; ",
+   "Wave 3: ", fmt(pc_w3_ctrl["mean"]), " [", fmt(pc_w3_ctrl["sd"]), "], N=", fmt0(pc_w3_ctrl["N"]), ". "
+ )
> 
> cat("\n[Figure 3 — First PC note]\n",     caption_pc,    "\n", sep = "")

[Figure 3 — First PC note]
Control (Entertainment) means [SD], N — Wave 1: -0.01 [0.95], N=391; Wave 2: -0.07 [0.97], N=348; Wave 3: -0.11 [0.95], N=313. 
> 
> 
> ###############################################
> # Appendix C - Persuasive effects by subgroups (media preference)
> ###############################################
> 
> w123$pid3w1_leaner <- ifelse(w123$party1 == 1, 'Democrat',
+                              ifelse(w123$party1 == 2, 'Republican',
+                                     ifelse(w123$party1 %in% c(3, 4) & w123$party4 == 1, 'Republican',
+                                            ifelse(w123$party1 %in% c(3, 4) & w123$party4 == 2, 'Democrat',
+                                                   ifelse(w123$party1 %in% c(3, 4) & w123$party4 == 3, 'Independent', NA)))))
> 
> w123$pid3w1_leaner <- factor(w123$pid3w1_leaner, levels = c('Democrat', 'Republican', 'Independent'))
> 
> 
> 
> # Compute the cross table
> cross_table <- table(w123$pid3w1_leaner, w123$med_choice)
> latex_table <- xtable(cross_table, 
+                       caption = "Cross Table of Party ID and Media Preference",
+                       label = "tab:cross_table")
> print(latex_table, type = "latex", include.rownames = TRUE)
% latex table generated in R 4.5.2 by xtable 1.8-4 package
% Fri Jan 30 18:36:36 2026
\begin{table}[ht]
\centering
\begin{tabular}{rrrr}
  \hline
 & Entertainment & Fox & MSNBC \\ 
  \hline
Democrat & 424 & 328 & 611 \\ 
  Republican & 351 & 773 & 152 \\ 
  Independent & 161 & 132 &  90 \\ 
   \hline
\end{tabular}
\caption{Cross Table of Party ID and Media Preference} 
\label{tab:cross_table}
\end{table}
> print(latex_table, type = "latex", file = "Output/tableA3_cross_table.tex", include.rownames = TRUE)
> 
> 
> # First Principal Component Plot
> pc_w1_fox <- lm(pcw1 ~ msnbc + fox, data = filter(w123,forcedchoice==1 & forcedchoice_w2==0 & med_pref == "Fox"))
> pc_w2_fox <- lm(pcw2 ~ msnbc + fox, data = filter(w123,forcedchoice==1 & forcedchoice_w2==0 & med_pref == "Fox"))
> pc_w3_fox <- lm(pcw3 ~ msnbc + fox, data = filter(w123,forcedchoice==1 & forcedchoice_w2==0 & med_pref == "Fox"))
> pc_w1_msnbc <- lm(pcw1 ~ msnbc + fox, data = filter(w123,forcedchoice==1 & forcedchoice_w2==0 & med_pref == "MSNBC"))
> pc_w2_msnbc <- lm(pcw2 ~ msnbc + fox, data = filter(w123,forcedchoice==1 & forcedchoice_w2==0 & med_pref == "MSNBC"))
> pc_w3_msnbc <- lm(pcw3 ~ msnbc + fox, data = filter(w123,forcedchoice==1 & forcedchoice_w2==0 & med_pref == "MSNBC"))
> pc_w1_ent <- lm(pcw1 ~ msnbc + fox, data = filter(w123,forcedchoice==1 & forcedchoice_w2==0  & med_pref == "Entertainment"))
> pc_w2_ent <- lm(pcw2 ~ msnbc + fox, data = filter(w123,forcedchoice==1 & forcedchoice_w2==0 & med_pref == "Entertainment"))
> pc_w3_ent <- lm(pcw3 ~ msnbc + fox, data = filter(w123,forcedchoice==1 & forcedchoice_w2==0 & med_pref == "Entertainment"))
> 
> marijuana_ests_fox <- tidy(pc_w1_fox) %>% mutate(Wave = 1,`Media Preference` = "Stated Preference: Fox") %>%
+   bind_rows(tidy(pc_w2_fox) %>% mutate(Wave = 2,`Media Preference` = "Stated Preference: Fox")) %>%
+   bind_rows(tidy(pc_w3_fox) %>% mutate(Wave = 3,`Media Preference` = "Stated Preference: Fox"))
> marijuana_ests_msnbc <- tidy(pc_w1_msnbc) %>% mutate(Wave = 1,`Media Preference` = "Stated Preference: MSNBC") %>%
+   bind_rows(tidy(pc_w2_msnbc) %>% mutate(Wave = 2,`Media Preference` = "Stated Preference: MSNBC")) %>%
+   bind_rows(tidy(pc_w3_msnbc) %>% mutate(Wave = 3,`Media Preference` = "Stated Preference: MSNBC"))
> marijuana_ests_ent <- tidy(pc_w1_ent) %>% mutate(Wave = 1,`Media Preference` = "Stated Preference: Entertainment") %>%
+   bind_rows(tidy(pc_w2_ent) %>% mutate(Wave = 2,`Media Preference` = "Stated Preference: Entertainment")) %>%
+   bind_rows(tidy(pc_w3_ent) %>% mutate(Wave = 3,`Media Preference` = "Stated Preference: Entertainment"))
> 
> marijuana_ests_all <- marijuana_ests_fox %>% 
+   bind_rows(marijuana_ests_msnbc) %>%
+   bind_rows(marijuana_ests_ent) %>%
+   mutate(`Media Preference` = factor(`Media Preference`,levels=c("Stated Preference: MSNBC","Stated Preference: Entertainment","Stated Preference: Fox"),ordered=T),
+          upper = estimate + qnorm(0.975)*std.error,
+          lower = estimate + qnorm(0.025)*std.error,
+          upper_90 = estimate + qnorm(0.95)*std.error,
+          lower_90 = estimate + qnorm(0.05)*std.error,
+          term = dplyr::recode(term,"msnbc"="MSNBC","fox" = "Fox"))
> 
> labs <- tibble(
+   Wave = c(1.5, 1.5), 
+   estimate = c(0.4, -0.3),
+   label = c("Fox rather\nthan Entertainment", "MSNBC rather\nthan Entertainment"),
+   term = c("Fox", "MSNBC"),
+   `Media Preference` = "Stated Preference: MSNBC"
+ ) %>%
+   mutate(`Media Preference` = factor(`Media Preference`,
+                                      levels = c("Stated Preference: MSNBC", "Stated Preference: Entertainment", "Stated Preference: Fox"),
+                                      ordered = T))
> 
> pdf('Output/figA1_pc.pdf', width=6, height=8)
> print(ggplot(filter(marijuana_ests_all,term!="(Intercept)"), aes(y=estimate, x=Wave, color=term)) +
+   geom_hline(yintercept = 0,lty=2) + 
+   geom_point(position=position_dodge(width=0.25))+ 
+   geom_line(aes(group=term),position=position_dodge(width=0.25)) + 
+   geom_errorbar(aes(ymin=lower, ymax=upper,width=0),size=0.5,position=position_dodge(width=0.25)) + 
+   geom_errorbar(aes(ymin=lower_90, ymax=upper_90,width=0),size=1,position=position_dodge(width=0.25)) + 
+   facet_wrap(~`Media Preference`,ncol=1) + 
+   geom_text(data=labs,aes(label=label)) +
+   theme_bw() + 
+   scale_x_continuous('Survey Wave',breaks=c(1,2,3),limits=c(0.75,3.25)) + 
+   scale_y_continuous('Effect of Treatment on\nFirst Principal Component',
+                      breaks=seq(-0.6,0.6,0.3),
+                      labels=c("\n\n-0.6\n(More\nliberal)", "-0.3",  "0.0","0.3", "(More\nconservative)\n0.6\n\n"),limits=c(-0.72,0.72)) + 
+   scale_color_manual("Treatment",breaks = c("Fox","MSNBC"),values=c(red_mit,blue_mit)) +
+   theme(legend.position="none",strip.background = element_rect(fill="lightgrey"),axis.title.y = element_text(margin = margin(r = 0))))
> dev.off()
null device 
          1 
> 
> 
> 
> ###############################################
> # Figure 6 - Accumulation effects
> ###############################################
> w123 <- w123 %>%
+   mutate(dosage = case_when(msnbc_w2==1 & msnbc==1 ~ "MSNBC twice",
+                             msnbc_w2==1 & entertainment==1 ~ "MSNBC once",
+                             fox_w2==1 & fox==1 ~ "Fox twice",
+                             fox_w2==1 & entertainment==1 ~ "Fox once",
+                             entertainment_w2==1 & entertainment==1 ~ "Entertainment twice"
+   ),
+   dosage = relevel(factor(dosage),ref = "Entertainment twice"))
> 
> 
> # First Principal Component Plot
> pc_w2_omnibus <- lm(pcw2 ~ dosage, data = filter(w123,forcedchoice==1 & forcedchoice_w2==1))
> 
> marijuana_ests <- tidy(pc_w2_omnibus) %>%
+   mutate(upper = estimate + qnorm(0.975)*std.error,
+          lower = estimate + qnorm(0.025)*std.error,
+          upper_90 = estimate + qnorm(0.95)*std.error,
+          lower_90 = estimate + qnorm(0.05)*std.error,
+          term = str_replace(term,"dosage",""),
+          dosage = str_replace(term,"Fox ",""),
+          dosage = str_replace(dosage,"MSNBC ",""),
+          term = str_replace(term," once",""),
+          term = str_replace(term," twice",""),
+          term = factor(term,levels=c("(Intercept)","Fox","MSNBC"),ordered=T)
+   )
> 
> pdf('Output/fig6_pc.pdf', width=6, height=4)
> print(ggplot(filter(marijuana_ests,term != "(Intercept)"), aes(y=estimate, color=term, alpha=dosage,shape=dosage, x=1)) +
+   geom_hline(yintercept = 0,lty=2) + 
+   geom_point(size=3,position=position_dodge(width=0.2)) + 
+   # geom_line(aes(group=term),position=position_dodge(width=0.25)) + 
+   geom_errorbar(aes(ymin=lower, ymax=upper,width=0),size=0.5,position=position_dodge(width=0.2)) + 
+   geom_errorbar(aes(ymin=lower_90, ymax=upper_90,width=0),size=1,position=position_dodge(width=0.2)) + 
+   facet_wrap(~term,scales = "free_x") + 
+   # geom_text(data=labs,aes(label=label)) +
+   theme_bw() + 
+   scale_x_continuous('Waves Treated with Partisan Media',breaks=c(0.95,1.05),labels=c("Wave 2\nonly","Wave 1\nand 2")) + 
+   scale_y_continuous('Effect of Treatment on\nFirst Principal Component',
+                      breaks=seq(-0.5,0.5,0.25),
+                      labels=c("\n\n-0.5\n(More\nliberal)", "-0.25",  "0.0",  "0.25",  "(More\nconservative)\n0.5\n\n"),limits=c(-0.51,0.51)) + 
+   scale_color_manual("Treatment",breaks = c("Fox","MSNBC"),values=c(red_mit,blue_mit)) +
+   scale_alpha_discrete("Wave 1 Treatment",range=c(0.25,3)) + 
+   scale_shape_manual("Wave 1 Treatment",values=c("once"=17,"twice"=16)) + 
+   theme(legend.position="none",strip.background = element_rect(fill="lightgrey"),axis.title.y = element_text(margin = margin(r = -20))))
Warning message:
Using alpha for a discrete variable is not advised. 
> dev.off()
null device 
          1 
> 
> 
> # Effect ratio estimates 
> 
> save_ratio_tex <- function(ratio, filename) {
+   formatted_ratio <- format(round(ratio, 2))
+   writeLines(formatted_ratio, paste0("Output/Estimates/", filename, ".tex"))
+ }
> coef_summary <- summary(pc_w2_omnibus)$coefficients
> 
> fox_once_effect <- coef_summary["dosageFox once", "Estimate"]
> fox_twice_effect <- coef_summary["dosageFox twice", "Estimate"]
> fox_twice_se <- coef_summary["dosageFox twice", "Std. Error"]
> fox_ratio_lower <- (fox_twice_effect - 1.96 * fox_twice_se) / fox_once_effect
> fox_ratio_upper <- (fox_twice_effect + 1.96 * fox_twice_se) / fox_once_effect
> save_ratio_tex(fox_ratio_lower, "fox_two_dose_ratio_lower")
> save_ratio_tex(fox_ratio_upper, "fox_two_dose_ratio_upper")
> 
> 
> msnbc_once_effect <- coef_summary["dosageMSNBC once", "Estimate"]
> msnbc_twice_effect <- coef_summary["dosageMSNBC twice", "Estimate"]
> msnbc_twice_se <- coef_summary["dosageMSNBC twice", "Std. Error"]
> msnbc_ratio_lower <- (msnbc_twice_effect - 1.96 * msnbc_twice_se) / msnbc_once_effect
> msnbc_ratio_upper <- (msnbc_twice_effect + 1.96 * msnbc_twice_se) / msnbc_once_effect
> save_ratio_tex(msnbc_ratio_lower, "msnbc_two_dose_ratio_lower")
> save_ratio_tex(msnbc_ratio_upper, "msnbc_two_dose_ratio_upper")
> 
> print(paste("Fox ratio range:", round(fox_ratio_lower, 3), "to", round(fox_ratio_upper, 3)))
[1] "Fox ratio range: -4.968 to 5.649"
> print(paste("MSNBC ratio range:", round(msnbc_ratio_lower, 3), "to", round(msnbc_ratio_upper, 3)))
[1] "MSNBC ratio range: 4.922 to -2.287"
> 
> 
> ###############################################
> # Figure 7 - Minimum Sample Size
> ###############################################
> 
> calc_mss <- function(effect, power=.8, attrition=0.12, v, alpha=.05, frac_used=1/9){
+   M <- qnorm(1-alpha/2, lower.tail=TRUE) + qnorm(power, lower.tail = TRUE)
+   (M *  sqrt(2*v/((1-attrition) * .25*(frac_used))))^2/effect^2
+ }
> 
> make_mss_df <- function(baseline_multiple, multiples, frac_used){
+   mss <- data.frame(multiples, 
+                     Fox = unlist(lapply(effects_fox*baseline_multiple, calc_mss, v=var(w123$pcw1, na.rm=TRUE), frac_used=frac_used)),
+                     Fox_lower = unlist(lapply(upper_fox*baseline_multiple, calc_mss, v=var(w123$pcw1, na.rm=TRUE), frac_used=frac_used)),
+                     MSNBC = unlist(lapply(effects_msnbc*baseline_multiple, calc_mss, v=var(w123$pcw1, na.rm=TRUE), frac_used=frac_used)),
+                     MSNBC_lower = unlist(lapply(upper_msnbc*baseline_multiple, calc_mss, v=var(w123$pcw1, na.rm=TRUE), frac_used=frac_used)))
+   
+   mss$multiples <- mss$multiples + 1
+   mss$baseline_multiple <- baseline_multiple
+   mss$frac_used <- frac_used
+   mss
+ }
> 
> make_mss_df_base_effect <- function(frac_used, baseline_multiples, multiples){
+   do.call(rbind, lapply(baseline_multiples, make_mss_df, multiples= multiples, frac_used=frac_used))
+ }
> 
> mod <- lm(pcw1 ~ fox + msnbc, data = filter(w123, forcedchoice == 1 & forcedchoice_w2 == 0))
> multiples <- seq(0, 3, by=.05)
> 
> effects_msnbc <- rep(coef(mod)[["msnbc"]], length(multiples)) * multiples  
> upper_msnbc <- rep(coef(mod)[["msnbc"]], length(multiples)) * multiples - 1.96 * sqrt(diag(vcov(mod))["msnbc"])
> 
> effects_fox <- rep(coef(mod)[["fox"]], length(multiples)) * multiples  
> upper_fox <- rep(coef(mod)[["fox"]], length(multiples)) * multiples + 1.96 * sqrt(diag(vcov(mod))["fox"])
> 
> fracs <- c(1/9, 1/5, 1/3)
> mss <- do.call(rbind, lapply(fracs, make_mss_df_base_effect, baseline_multiples=1:2, multiples=multiples))
> mss$n_groups <- paste0(mss$frac_used^(-1), ' Treatment Groups')
> 
> names(mss)[names(mss) == "multiples"] <- "Gamma"
> 
> mss_long <- pivot_longer(mss[, c('Gamma', 'Fox', "MSNBC", "baseline_multiple", "n_groups")], 
+                          cols=c('Fox','MSNBC'), 
+                          names_to = "name", 
+                          values_to = "value")
> mss_long$baseline_multiple <- ifelse(mss_long$baseline_multiple == 1, 'Empirical Effect', '2X Empirical Effect')
> 
> pdf('Output/fig7_mss.pdf', width=10, height=5)
> print(ggplot(mss_long[mss_long$Gamma <= 2.5 & mss_long$Gamma > 1.50 & mss_long$baseline_multiple=='Empirical Effect',], 
+        aes(x=Gamma, y=value, color=n_groups)) + 
+   geom_point() + 
+   geom_line() + 
+   geom_hline(yintercept=4886, alpha=.5) + 
+   geom_hline(yintercept=0, alpha=.25, linetype=2) +
+   facet_grid(~name) +
+   theme_bw() + 
+   theme(legend.title = element_blank(),
+         strip.background = element_rect(fill="lightgrey"), 
+         axis.title.y = element_text(margin = margin(r = 20)),
+         plot.background = element_rect(fill = 'white'), 
+         panel.background = element_rect(fill = 'white')) + 
+   scale_x_continuous(labels = scales::percent) +
+   scale_y_continuous(labels = scales::comma) +
+   labs(x = 'Effect of Two Doses as % of Effect of One',
+        y = 'Minimum Sample Size (MSS)'))
> dev.off()
null device 
          1 
> 
> 
> ###############################################
> # Appendix B - Power Calculation
> ###############################################
> 
> table_data <- mss %>%
+   filter(baseline_multiple == 1, 
+          n_groups == "9 Treatment Groups", 
+          Gamma %in% seq(1.25, 4, by = 0.25))
> table_data <- table_data[, c("Gamma", "Fox", "Fox_lower", "MSNBC", "MSNBC_lower")]
> 
> table_data <- table_data %>%
+   mutate(across(c(Fox, Fox_lower, MSNBC, MSNBC_lower), round))
> 
> table_a1 <- kable(table_data, 
+                   format = "latex",
+                   col.names = c("$\\gamma$", "MSS", "95\\% CI Lower Bound", "MSS", "95\\% CI Lower Bound"),
+                   caption = "Minimum Sample Size (MSS) Needed for 80\\% Power By Proportional Change in Treatment Effect ($\\gamma$) Based on Wave 1 Treatment",
+                   align = c('c', 'r', 'r', 'r', 'r'),
+                   booktabs = TRUE,
+                   escape = FALSE) %>%
+   kable_styling(latex_options = c("striped", "hold_position")) %>%
+   add_header_above(c(" " = 1, "Fox" = 2, "MSNBC" = 2))
> 
> writeLines(table_a1, "Output/tableA2.tex")
> 
> ###############################################
> # Save estimates
> ###############################################
> 
> save_estimate_tex <- function(estimate, filename) {
+   formatted_estimate <- format(round(estimate, 3))
+   writeLines(formatted_estimate, paste0("Output/Estimates/", filename, ".tex"))
+ }
> 
> extract_and_save_estimates <- function(model, wave, outcome) {
+   coef_data <- coef(summary(model))
+   treatments <- c("fox", "msnbc")
+   
+   for (treatment in treatments) {
+     estimate <- coef_data[treatment, "Estimate"]
+     std_error <- coef_data[treatment, "Std. Error"]
+     save_estimate_tex(estimate, paste0("wave", wave, "_", outcome, "_", treatment, "_effect"))
+     save_estimate_tex(estimate - 1.96 * std_error, paste0("wave", wave, "_", outcome, "_", treatment, "_effect_lower"))
+     save_estimate_tex(estimate + 1.96 * std_error, paste0("wave", wave, "_", outcome, "_", treatment, "_effect_upper"))
+   }
+ }
> 
> extract_and_save_estimates(pc_w1, 1, "marijuana_pc")
> extract_and_save_estimates(pc_w2, 2, "marijuana_pc")
> extract_and_save_estimates(pc_w3, 3, "marijuana_pc")
> 
> proc.time()
   user  system elapsed 
  1.630   0.105   1.697 
